The Breaking Strain of Neutron Star Crust and Gravitational Waves 
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Mountains on rapidly rotating neutron stars efficiently radiate gravitational waves. The maximum 
possible size of these mountains depends on the breaking strain of neutron star crust. With multi- 
million ion molecular dynamics simulations of Coulomb solids representing the crust, we show that 
the breaking strain of pure single crystals is very large and that impurities, defects, and grain 
boundaries only modestly reduce the breaking strain to around 0.1. Due to the collective behavior 
of the ions during failure found in our simulations, the neutron star crust is likely very strong 
and can support mountains large enough so that their gravitational wave radiation could limit the 
spin periods of some stars and might be detectable in large scale interferometers. Furthermore, 
our microscopic modeling of neutron star crust material can help analyze mechanisms relevant in 
magnetar giant and micro flares. 



"Mountains" on rapidly rotating neutron stars (NS) 
efficiently radiate gravitational waves (GW) because of 
their large masses and high accelerations[l|. Searches 
for these waves are presently underway with large scale 
interferometers 0, Q. Furthermore, the angular momen- 
tum radiated could control the spin periods of some ac- 
creting NS^]. As material falls on a NS, from its binary 
companion, the angular momentum accreted spins up the 
star. However, GW radiation increases very rapidly with 
rotational frequency. The star may reach spin equilib- 
rium where the angular momentum gained from accretion 
is balanced by that lost to GW radiation. This could ex- 
plain why the most rapidly spinning pulsars observed are 
only spinning at about half of the breakup rate. 

How high can a mountain be, before it collapses under 
the NS's extreme gravity? This depends on the breaking 
strain (BS)[i] of NS crust (NSC) which is the fractional 
deformation when the crust fails. Estimates of the break- 
ing strain vary by orders of magnitude and are based on 
lose analogies with conventional materials rather than 
on detailed calculations or simulations^, This is, per- 
haps, the largest uncertainty in estimating upper limits 
for these GW. In this paper we use molecular dynamics 
(MD) simulations to significantly improve estimates of 
the breaking strain. 

Crust breaking may also be important for magnetar 
giant flares. These are extremely energetic Gamma-ray 
bursts from strongly magnetized NS .8J. Thompson and 
Duncan model giant flares as "star quakes" @ . A strong 
twisted magnetic field stresses the crust until it breaks. 



The crust then moves and allows magnetic field lines to 
reconnect. This releases magnetic energy that powers 
the flare. This mechanism may require the crust to be 
relatively strong in order for it to control the very strong 
magnetic field. 

How the crust breaks may be important for the ex- 
citation of NS oscillation modes from a "star quake". 
Quasiperiodic oscillations have been observed in the tails 
of giant flares EH- These have been interpreted as 
shear oscillations of the crust. If these modes can be 
convincingly identified, they may provide considerable 
information on the structure of the NS and its crust fl^. 
In addition, LIGO (Laser Interferometer Gravitational 
Wave Observatory) has searched for GW from a mag- 
netar giant fiare that could come from large amplitude 
NS oscillations[13]. Our simulations of crust breaking 
may lead to insight into the amplitudes that might be 
expected. 

In the crust of a NS, electrons form a very degenerate 
relativistic gas. The ions are completely pressure ion- 
ized and have Coulomb interactions that are screened at 
large distances by the slightly polarizable electron gas. 
The interaction between two ions i,j is assumed to be a 
screened Coulomb or Yukawa potential 



2^2 



-rij/\a 



(1) 



'Electronic address; lhorowit@indiana.edul 
t Electronic address: | kkadau@lanl.gov] 



where the ions have charge Z, Vij is the distance be- 
tween them, and the electron screening length Ag = 
7r^/^/2e(37r^ne)^/'^ with Tie the electron density. The to- 
tal potential energy is given by the sum over all pairs 
Tlii^j i'irij)- Charge neutrality insures that = Zn 
where n is the ion density. The ions are assumed to form 
a classical one component plasma that can be character- 
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ized by the Coulomb parameter F, 



r = 
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This is the ratio of a typical Coulomb energy to thermal 
energy and the ion sphere radius a = [3/(47™)]^/"^ is a 
typical distance between ions. 

We calculate at a reference density ofn = 7.18xl0~^ 
fm^'' using Z = 29.4 and an atomic mass number A — 
88 (1 X 10^'^ g/cm'^)[fH|. Alternatively, nuclei are very 
neutron rich near the base of the crust. For A = 880 this 
n corresponds to 10^** g/cm'^. Results can be scaled to 
other densities at constant F and approximately scaled to 
other Z at constant F, which involves only a small change 
in Ag. Most of our simulations are for a temperature 
T = 0.1 MeV. This corresponds to F = 834. We ignore 
any strong interactions between ions and the effects of 
free neutrons that are present in the inner crust. 
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FIG. 1: Shear stress versus strain for perfect and defective 
bcc and fee single crystals containig about 2 million ions for 
different shear systems — as given by the crystallographic ori- 
entation of the shear plane and the direction of the applied 
stress — are shown. In addition a poly crystalline sample 
containing 12.8 million ions and 8 randomly oriented grains 
with an average grain diameter of 3962 fm is shown. Results 
were obtained at a strain rate of 4x10^^ c/fm. 

We perform large-scale MD simulations of shearing 
NSC material as we believe that this is the most impor- 
tant mode of failure in the crust, although real neutron 
star strain fields could involve some combination of shear 
strain and tension. Tension simulations at constant vol- 
ume, had a breaking strain and strength that was smaller 
by 2.5 and 2, respectively. We use two independent MD 
codes: YukawaMD is a serial code where the system is 
strained by deforming the periodic boundary. An origi- 
nally cubic simulation volume is deformed according to 
X ^ X + ey/2, y ^ y + ex/2, and z {1 + e'^/2)z. Here 
the strain e is increased linearly with time. 



SPaSMfiej is a high-performance parallel MD code 
where the system is strained by moving a boundary layer 
of frozen ions on top and bottom against each other [171. 
In order to simulate systems containing up to 12.8 mil- 
lion ions we employed up to 512 processors of a paral- 
lel architecture consisting of AMD Opteron cores @ 2.2 
GHz installed at Los Alamos ( Yellowrail) . The computa- 
tional cost for these simulations are about 100 times more 
as compared to short range Lennard-Jones computations 
because of the much longer ranged force (rc = lOAe). 
The shear is initiated by having a frozen layer of ions, 
i.e. ions that are not dynamic, on the top and bottom 
of the system. The layers are separated by the distance 
d and moved with a velocity v against each other to ini- 
tiate the desired shear rate de/dt = v/d. We cutoff the 
potential at the large distance Vc — 10Ae|18j resulting 
in each ion interacting with about 5400 neighbors; for 
smaller values of Tc the results depend on rc[l^. See 
supplemental Figures (S1,S2) 

Figure [T] shows the shear stress versus strain for body- 
centered cubic (bcc) crystals, as this is the equilibrium 
structure. Some results for the face-centered cubic (fee) 
structures are added for comparison. For all investigated 
crystallographic shear systems the perfect crystals show 
a BS well above 0.1, and break in a rather abrupt fash- 
ion with only a very small region where plasticity, i.e. 
deviation from a linear stress-strain relation, is present. 
The multi-million ion systems were strained at a rate 
of 4x10"'' c/fm. As we cannot simulate the very large 
time and length scales associated with NSC we have to 
rely on estimates based on a series of simulations that 
suggested basically no size effects for the single crystal 
simulations presented here (Fig. S3 [20]) and a converg- 
ing result at low shear rates (Fig. S4 |20j). We also note 
that tripling the temperature only reduces the maximum 
stress by about 25% and does not significantly alter the 
BS. 

As very little is known about the defect structure, grain 
sizes, and associated grain boundaries of NSC that po- 
tentially can reduce the strength, we also consider poly- 
crystalline materials that are generated using a Voronoi 
construction in which grain centers and orientations are 
picked at random and each ion belongs to the grain center 
that is closest[2l|. After this initialization, the system is 
equilibrated by heating it from 0.1 to 0.3 MeV and back 
to 0.1 MeV over a total simulation time of 275000 fm/c. 
During equilibration at elevated temperatures, ions find 
lower energy configurations exponentially faster, relax- 
ing the initial artificial grain boundary structures more 
rapidly. The system is then sheared at a strain rate of 
4 X 10^^ c/fm. Figure 2 shows one system with 12.8 mil- 
lion ions at strains of 0.0, 0.05, 0.1, and 0.15 (for stress- 
strain curve see Fig.l). The radial distribution func- 
tions exhibit characteristic peaks of the bcc structure for 
strains up to 0.1. After failure, at 0.15 strain some degree 
of amorphization might be possible. The system starts 
to deform plastically near the grain boundaries in a col- 
lective manner without exhibiting signs of dislocations or 
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FIG. 2: Poly crystalline sample with 12.8 million ions con- 
sisting of 8 differently oriented grains with an average grain 
diameter of 3961 fm at strains of 0.0, 0.05, 0.1, and 0.15. Up- 
per panel: The radial distribution functions at different shear 
states exhibit the characteristic peaks of the bcc structure for 
strains up to 0.1. At 0.15 strain some degree of amorphization 
might be possible. Lower panel: At 0.0 strain the eight dif- 
ferent grains are shown in different colors. At strains of 0.05, 
0.1, and 0.15 the red color indicates plastic deformation, i.e. 
deviations of the ions away from an ideally uniformly sheared 
bcc lattice in the range of to two times the nearest neighbor 
distance. Note that the top and bottom layers are frozen and 
moved to impose the shear. 



other more localized events (Fig. 2, See also supplemental 
Movies S5,S6 However, at large strains the samples 

are plastically deformed and the mode of failure seems to 
be a reorientation of local regions — possibly associated 
with some degree of amorphization (Fig. 2) — to accom- 
modate the shear. This collective, rather than more lo- 
calized dislocation-based, mode of plasticity makes the 
crystal stronger and break at larger strains than terres- 
trial metals, where the electronic density can have a lo- 



calized structure to accommodate local defects. This fail- 
ure mechanism also does not allow for voids or fracture 
to appear as these localized defects would possibly heal 
under the influence of the high pressure. In fact, sim- 
ulations that started out with a cylindrical hole with a 
diameter of 2.5 x the nearest neighbor distance initial- 
ized into the otherwise perfect single crystal, showed the 
void quickly heals under the influence of the enormous 
pressure present in the system and the maximum stress 
is only reduced by 25% (Fig.l, Movie S7 [13]). This is 
consistent with the prediction of Jones that voids would 
not form because of the high pressure [2^. Furthermore, 
simulations containing only two grains and no frozen-in 
ions rapidly evolved into a single crystal [l^. We investi- 
gated samples with grain diameters between 570 fm and 
3961 fm that showed a clear monotonic trend to increase 
the BS and stress with the grain size. As we expect to 
have larger grains in NSC material than we could possi- 
bly simulate on a computer we expect the crust material 
strength only to be slightly reduced by the presence of 
grain boundaries. Also, note that accreting NS can form 
new crust very slowly over accretion times of thousands 
of years. This should allow plenty of time for larger crys- 
tals to form. Second, samples with very small grain size 
may have a reduced thermal conductivity and this could 
disagree with crust cooling observations [23, ^25!]. Indeed, 
the grain size could be much larger than in our simula- 
tions. In that case, we expect the strength of the system 
to be larger and closer to our single crystal results. 

To study the role of impurities we form a 27648 ion 
crystal made of the complex rapid proton capture nucle- 
osynthesis ash considered in ref.flSl. We use the solid 
composition from Table I of ref.[15]. This includes a 
range of ions from Z = 8 to 47 and can be character- 
ized by an impurity parameter Q = (Z^) — {Z)"^ = 22.3 
that describes the dispersion in ion charges. Shternin 
et al.[2^ and Brown and Gumming [23| study the rapid 
crust cooling of two NS, after extended outbursts. Brown 
and Gumming conclude that Q must be less than 10 24| . 
see also[23]. Otherwise impurities would reduce the crust 
thermal conductivity by too much to agree with obser- 
vation. Our crystal has more impurities than this limit. 
Therefore, we believe our results may represent an up- 
per limit on the effects of impurities. Indeed, we have 
so many impurities that our runs with significantly more 
impurities phase separate 15[ . We equilibrate our crystal 
by heating it up and then slowly cooling it back down to 
T = 0.1 MeV over a total simulation time of about 85 
million fm/c. We shear this crystal at a strain rate of 
8xl0"^c/fm and find a maximum stress that is reduced 
from the maximum stress of a pure crystal by less than 
45% at a similar strain. We conclude that impurities are 
unlikely to significantly reduce the strength of NSG. 

Our simulations show that the breaking strength of 
neutron star crust is about 10^° times more than for ter- 
restrial engineering materials such as metal-alloys where 
the strength is measured in fractions of a GPa. The 
largest contributor to this tremendous difference is of 



4 



course the enormous pressure and thus density of the 
crust. Furthermore, the screened Coulomb interaction is 
purely repulsive and has no explicit length scalefl^, i.e. 
the system at twice the density behaves just like the sys- 
tem at the original density only at a lower temperature 
(Eql2|). This causes the material to fail abruptly in a 
collective manner at a large strain, rather than yielding 
continuously at low strain as observed in metals, because 
of the formation of dislocations. For example, the break- 
ing strain of steel is around 0.005, some twenty times 
smaller than what we find for the neutron star crust. We 
speculate that the collective plastic behavior found here 
could help to improve design strategies that suppress the 
weakening effects of dislocations and other more local- 
ized defects in conventional materials. Note that small 
Coulomb solids have been studied in the laboratory using 
cold trapped ions[27j. 

Gravitational wave radiation depends on the elliptic- 
ity e (fractional difference in moments of inertia) of a 
rotating star. We estimate, following ref. that our 
breaking strain w 0.1 can support an e < 4 x 10~^ for 
a 1.4 solar mass, 10 km radius NS. This ellipticity, for 
a rapidly rotating star, will generate GW that could be 
detectable by LIGO Q. 

In conclusion, we have performed large-scale MD sim- 
ulations of the breaking strain (BS) of Coulomb solids 



representing neutron star crust (NSC). We find a collec- 
tive mode of failure that does not involve localized defects 
such as dislocations, opening of voids, or fractures. The 
breaking strain in the presence of introduced defects and 
impurities is only moderately reduced to about 0.1. The 
large breaking strain, that we find, should support moun- 
tains on rapidly rotating NS large enough to efhciently 
radiate gravitational waves. This motivates further work 
on mountain forming mechanisms and searches for con- 
tinuous gravitational waves. Furthermore, the methods 
we presented here are very promising to describe other 
aspects of Coulomb solids and neutron star crust in par- 
ticular. 
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